#include <string>
#include <vector>
#include <cstdlib>
#include <iomanip>
#include <fstream>
#include <iostream>
#include <assert.h>
#include <stdlib.h>
#ifdef __USE_PTHREADS__
#include <pthread.h>
#endif
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#ifdef __CPUTIME__
#define CLKS_T std::vector< timespec * > &
#else
#define CLKS_T timespec *
#endif
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#include <Types.h>
#include <Kernels.h>
#include <Defines.h>
#include <TimeDiff.h>
#include <DagLoader.h>
#include <QTSKernel.h>
#include <MallocAlign.h>
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#include <errno.h>
#define handle_error_en(en, msg) do { errno = en; perror(msg); exit(EXIT_FAILURE); } while (0)
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
using namespace std;
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#ifdef __HUGE_PAGES__
#include <sys/shm.h>
int shmidNodes;
void HugeFree( int shmid ) {
   struct shmid_ds unused_buffer;
   printf( "Deleting shared segment" );
   int code = shmctl( shmid, IPC_RMID, &unused_buffer );
   printf( " -- return code %i\n", code );
}
#endif
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#ifdef LOCAL_DEBUG_ENABLED
#define __SERIALIZE_PTHREADS__
#else
// #define __SERIALIZE_PTHREADS__
#endif
#ifdef __USE_PTHREADS__
// #define __HL1__
#define __HL2__
// #define __HL4__
#ifdef __HL2__
const U32 numParts = 64;
const U32 numPartsPerKernel = 2;
#endif
#ifdef __HL4__
const U32 numParts = 256;
const U32 numPartsPerKernel = 4;
#endif
#else
const U32 numParts = 64;
#endif
const S32 numRoots   = numParts ==   4 ? 1 :
                       numParts ==  64 ? 2 :
                       numParts == 256 ? 4 : -1;
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
class COOArgs {
public:
   U32 numKernels;
   U32 tidx;
   U64 nnz;
   U64 m;
   U64 hgt;
   U64 * _nodes;
   U32 nnzValsRootNodePtr[ numRoots ];
   U32 colIdxsRootNodePtr[ numRoots ];
   const double * _x;
   double * _y;
   timespec * clks;
   bool hasLastKernel;
   
   COOArgs( void ) {
      numKernels = 0;
      hasLastKernel = false;
   }
   COOArgs( const COOArgs & o ) {
      numKernels = o.numKernels;
      tidx = o.tidx;
      nnz  = o.nnz;
      m    = o.m;
      hgt  = o.hgt;
      _nodes = o._nodes;
      _x = o._x;
      _y = o._y;
      clks = o.clks;
      for( S32 i = 0; i < numRoots; i++ ) {
         this->nnzValsRootNodePtr[ i ] = o.nnzValsRootNodePtr[ i ];
         this->colIdxsRootNodePtr[ i ] = o.colIdxsRootNodePtr[ i ];
      }
   }
   void operator =( const COOArgs & o ) {
      numKernels = o.numKernels;
      tidx = o.tidx;
      nnz  = o.nnz;
      m    = o.m;
      hgt  = o.hgt;
      _nodes = o._nodes;
      _x = o._x;
      _y = o._y;
      clks = o.clks;
      for( S32 i = 0; i < numRoots; i++ ) {
         this->nnzValsRootNodePtr[ i ] = o.nnzValsRootNodePtr[ i ];
         this->colIdxsRootNodePtr[ i ] = o.colIdxsRootNodePtr[ i ];
      }
   }
   void print( void ) const {
      cout << "COOArgs contents:" << endl;
      cout << "... numKernelLaunches = " << numKernels << endl;
      cout << "...  tidx = " << tidx << endl;
      cout << "...   nnz = " << nnz  << endl;
      cout << "...     m = " << m    << endl;
      cout << "...   hgt = " << hgt  << endl;
      cout << "... nodes = " << hex << _nodes << endl;
      cout << "... nnzValsRootNodePtr[ 0 ] = " << hex << nnzValsRootNodePtr[ 0 ] << endl;
      cout << "... colIdxsRootNodePtr[ 0 ] = " << hex << colIdxsRootNodePtr[ 1 ] << endl;
      cout << "...     x = " << hex << _x << endl;
      cout << "...     y = " << hex << _y << endl;
      cout << "...   clks = " << hex << clks << endl;
      cout << dec;
   }
};
struct DAGLevelPointers {

   U64 * nodes;
   S32 level0[   4 ];
   S32 level1[  16 ];
   S32 level2[  64 ];
   S32 level3[ 256 ];
   
};
struct MtxMetaData {
   U64 nnzHgt;
   U64 idxHgt;
   U32 m;
   U32 n;
   U32 nnz;
   DAGLevelPointers nnzRoots;
   DAGLevelPointers idxRoots;
};
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void PrintPTRVectorToMMFile( const double * vec, const U64 nnz, const bool complex, const std::string mtxFileName ) {
   
   U64 i;
   
   std::string dataType = "real";
   std::string symmType = "general";
   FILE * OFH;
   OFH = fopen( mtxFileName.c_str(), "w" );
   
   if( complex ) { dataType = "complex";   }

   #ifdef __MAC_OSX__
   fprintf( OFH, "%%MatrixMarket matrix coordinate %s %s\n", dataType.c_str(), symmType.c_str() );
   fprintf( OFH, "%llu 1 %llu\n", nnz, nnz );
   #else
   fprintf( OFH, "%%MatrixMarket matrix coordinate %s %s\n", dataType.c_str(), symmType.c_str() );
   fprintf( OFH, "%lu 1 %lu\n", nnz, nnz );
   #endif
   
   if( complex ) {
      for( i=0; i < nnz; i++ ) {
         U64 idx = i << 1ULL;
         double rv = vec[ idx + 0 ];
         double iv = vec[ idx + 1 ];
         #ifdef __MAC_OSX__
         fprintf( OFH, "%lld 1 %.25lf %.25lf\n", i+1, rv, iv );
         #else
         fprintf( OFH, "%ld 1 %.25lf %.25lf\n", i+1, rv, iv );
         #endif
      }
   }
   else {
      for( i=0; i < nnz; i++ ) {
         #ifdef __MAC_OSX__
         fprintf( OFH, "%lld 1 %.25lf\n", i+1, vec[i] );
         #else
         fprintf( OFH, "%ld 1 %.25lf\n", i+1, vec[i] );
         #endif
      }
   }
   fclose( OFH );
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
double GenRandDouble( void ) {
   long long int n;
   int    rn;
   double rf;
   double value;

   rn = rand();
   rf = ((double)rn) / ((double)RAND_MAX);
   rf *= 2e10;
   n = (long long int)rf;
   n -= (long long int)1e10;
   value = n/1e10;
   return value;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
inline void SetPThreadToRT_FIFO( void ) {

   int s;
   struct sched_param param;
   param.sched_priority = 80;

   s = pthread_setschedparam( pthread_self(), SCHED_FIFO, &param );
   if( s != 0 ) {
      handle_error_en(s, "pthread_setschedparam");
   }
}
void SetThreadAffinity( const U32 tidx ) {

   #ifndef __MAC_OSX__
   #ifdef __USE_PTHREADS__
   // notes for running on hcs001:
   // ... socket 0: 0, 1, 2, 3,  8,  9, 10, 11
   // ... socket 1: 4, 5, 6, 7, 12, 13, 14, 15
   // ... with centos:
   // ... ... prefers to run 8 threads distributed across both sockets
   // ... ... 10% performance boost from numactl --interleave all
   // ... ... can also use numactl -m [0|1]
   
   const U32 NUM_LOGICAL_CORES = 24;
   // const U32 cpuIdxsBySocket[ 16 ] = { 0,  8, 1,  9, 2, 10, 3, 11, 4, 12, 5, 13, 6, 14, 7, 15 };
   // const U32 cpuIdxsBySocket[ 16 ] = { 0,  1, 2,  3, 8, 9, 10, 11, 4, 5, 6, 7, 12, 13, 14, 15 };
   // const U32 cpuIdxsBySocket[ 16 ] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 };
   
   int s;
   int cpuIdx;
   cpu_set_t cpuset;
   pthread_t thread;
   
   thread = pthread_self();
   cpuIdx = tidx % NUM_LOGICAL_CORES;
   // cpuIdx = cpuIdxsBySocket[ cpuIdx ];
   
   // -- requires numa.h -- for( U32 i = 0; i < 16; i++ ) {
   // -- requires numa.h --    U32 socket = numa_node_of_cpu( i );
   // -- requires numa.h --    DEBUG_PRINT( "cpu: " << SETDEC( 2 ) << i << " is on socket " << SETDEC( 2 ) << socket );
   // -- requires numa.h -- }
   
   DEBUG_PRINT( "tidx = " << tidx << ", cpuIdx = " << cpuIdx );
   
   // set affinity mask
   CPU_ZERO( &cpuset );
   CPU_SET( cpuIdx, &cpuset );

   // set affinity
   s = pthread_setaffinity_np( thread, sizeof(cpu_set_t), &cpuset );
   if( s != 0 ) {
      handle_error_en( s, "pthread_setaffinity_np" );
      cout << endl;
      assert( false );
   }
   #endif
   #endif
}
inline void SetProcessSocketAffinity( U32 socketIdx ) {

   #ifndef __MAC_OSX__
   #ifdef __USE_PTHREADS__
   // Socket 0: ( 0 8 1 9 2 10 3 11 )
   // Socket 1: ( 4 12 5 13 6 14 7 15 )

   int s;
   cpu_set_t cpuset;
   // pid_t pid = getpid();
   
   // set affinity mask
   const U32 socket0[ 12 ] = { 0, 1, 2, 3,  4,  5, 12, 13, 14, 15, 16, 17 };
   const U32 socket1[ 12 ] = { 6, 7, 8, 9, 10, 11, 18, 19, 20, 21, 22, 23 };

   socketIdx = socketIdx % 2;
   const U32 * cpuIdxs = NULL;
   if( socketIdx == 0 ) { cpuIdxs = socket0; }
   if( socketIdx == 1 ) { cpuIdxs = socket1; }

   DEBUG_PRINT( "SetProcessSocketAffinity(): assigning process to socket " << socketIdx );   
   CPU_ZERO( &cpuset );
   for( U32 i = 0; i < 12; i++ ) {
      DEBUG_PRINT( "SetProcessSocketAffinity(): adding cpuIdx " << SETDEC( 2 ) << cpuIdxs[i] << " to cpuset mask" );   
      CPU_SET( cpuIdxs[i], &cpuset );
   }
   
   // set affinity
   // s = sched_setaffinity( pid, sizeof( cpu_set_t ), &cpuset );
   s = sched_setaffinity( 0, sizeof( cpu_set_t ), &cpuset );
   if( s != 0 ) {
      handle_error_en( s, "pthread_setaffinity_np" );
      cout << endl;
      assert( false );
   }
   #endif
   #endif
}
inline void SetSchedulerToRT_FIFO( void ) {

   #ifndef __MAC_OSX__
   int i;
   
   struct sched_param mysched;
   mysched.sched_priority = 80;  
   sched_setscheduler( 0, SCHED_FIFO, &mysched );
   
   i = sched_getscheduler(0);
   printf( "\n Scheduler %s\n", i == SCHED_FIFO ? "SCHED_FIFO" : "UNKNOWN" );
   #endif
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
inline void * SpMVFlatKernelCOO_HL4_last( void * vpargs ) {
   #ifdef __HL4__
   DEBUG_PRINT( "... ... entering SpMVFlatKernelCOO_HL4_last() " );
   COOArgs * pargs = (COOArgs*) vpargs;
   const U64 hgt = pargs->hgt;
   const U64 * _nodes = pargs->_nodes;

   const U32 nnzValsRootNodePtr0 = pargs->nnzValsRootNodePtr[ 0 ];
   const U32 colIdxsRootNodePtr0 = pargs->colIdxsRootNodePtr[ 0 ];
   const U32 nnzValsRootNodePtr1 = pargs->nnzValsRootNodePtr[ 1 ];
   const U32 colIdxsRootNodePtr1 = pargs->colIdxsRootNodePtr[ 1 ];
   const U32 nnzValsRootNodePtr2 = pargs->nnzValsRootNodePtr[ 2 ];
   const U32 colIdxsRootNodePtr2 = pargs->colIdxsRootNodePtr[ 2 ];
   const U32 nnzValsRootNodePtr3 = pargs->nnzValsRootNodePtr[ 3 ];
   const U32 colIdxsRootNodePtr3 = pargs->colIdxsRootNodePtr[ 3 ];
   
   const double * _x = pargs->_x;
   double * _y = pargs->_y;
   double * y0 = _y;
   double * y1 = _y;
   double * y2 = _y;
   double * y3 = _y;
   y0--;
   y1--;
   y2--;
   y3--;
   
   const U32 maxHgt   = 20;
   const U32 numAlloc = 4 * maxHgt;
   U32 lvl = hgt;

   U32 levelsStackX[ numAlloc ];
   U32 nnzPtrStack0[ numAlloc ];
   U32 idxPtrStack0[ numAlloc ];
   U32 nnzPtrStack1[ numAlloc ];
   U32 idxPtrStack1[ numAlloc ];
   U32 nnzPtrStack2[ numAlloc ];
   U32 idxPtrStack2[ numAlloc ];
   U32 nnzPtrStack3[ numAlloc ];
   U32 idxPtrStack3[ numAlloc ];
   
   const S32 * nodes = (const S32 *) _nodes;
   register S32 colIdx01 = 0;
   register S32 colIdx11 = 0;
   register S32 colIdx21 = 0;
   register S32 colIdx31 = 0;
   register double result0 = 0.0;
   register double result1 = 0.0;
   register double result2 = 0.0;
   register double result3 = 0.0;
   
   PLS upidx0;
   PLS upnnz0;
   PLS upidx1;
   PLS upnnz1;
   PLS upidx2;
   PLS upnnz2;
   PLS upidx3;
   PLS upnnz3;
   
   U32 * plvlX = &levelsStackX[0];
   upidx0.ps  = &idxPtrStack0[0];
   upnnz0.ps  = &nnzPtrStack0[0];
   upidx1.ps  = &idxPtrStack1[0];
   upnnz1.ps  = &nnzPtrStack1[0];
   upidx2.ps  = &idxPtrStack2[0];
   upnnz2.ps  = &nnzPtrStack2[0];
   upidx3.ps  = &idxPtrStack3[0];
   upnnz3.ps  = &nnzPtrStack3[0];
   const U32 * final = plvlX;
   
   const LSU2X * xnnzLSU0 = (const LSU2X *)&nodes[ nnzValsRootNodePtr0 ];
   const LSU2X * xidxLSU0 = (const LSU2X *)&nodes[ colIdxsRootNodePtr0 ];
   const LSU2X * xnnzLSU1 = (const LSU2X *)&nodes[ nnzValsRootNodePtr1 ];
   const LSU2X * xidxLSU1 = (const LSU2X *)&nodes[ colIdxsRootNodePtr1 ];
   const LSU2X * xnnzLSU2 = (const LSU2X *)&nodes[ nnzValsRootNodePtr2 ];
   const LSU2X * xidxLSU2 = (const LSU2X *)&nodes[ colIdxsRootNodePtr2 ];
   const LSU2X * xnnzLSU3 = (const LSU2X *)&nodes[ nnzValsRootNodePtr3 ];
   const LSU2X * xidxLSU3 = (const LSU2X *)&nodes[ colIdxsRootNodePtr3 ];

   *plvlX = lvl; plvlX++;
   *plvlX = lvl; plvlX++;
   *plvlX = lvl; plvlX++;
   *plvlX = lvl; plvlX++;
   
   *(upnnz0.pl) = xnnzLSU0->r.l; upnnz0.pl++;
   *(upnnz0.pl) = xnnzLSU0->l.l; upnnz0.pl++;
   *(upidx0.pl) = xidxLSU0->r.l; upidx0.pl++;
   *(upidx0.pl) = xidxLSU0->l.l; upidx0.pl++;

   *(upnnz1.pl) = xnnzLSU1->r.l; upnnz1.pl++;
   *(upnnz1.pl) = xnnzLSU1->l.l; upnnz1.pl++;
   *(upidx1.pl) = xidxLSU1->r.l; upidx1.pl++;
   *(upidx1.pl) = xidxLSU1->l.l; upidx1.pl++;

   *(upnnz2.pl) = xnnzLSU2->r.l; upnnz2.pl++;
   *(upnnz2.pl) = xnnzLSU2->l.l; upnnz2.pl++;
   *(upidx2.pl) = xidxLSU2->r.l; upidx2.pl++;
   *(upidx2.pl) = xidxLSU2->l.l; upidx2.pl++;

   *(upnnz3.pl) = xnnzLSU3->r.l; upnnz3.pl++;
   *(upnnz3.pl) = xnnzLSU3->l.l; upnnz3.pl++;
   *(upidx3.pl) = xidxLSU3->r.l; upidx3.pl++;
   *(upidx3.pl) = xidxLSU3->l.l; upidx3.pl++;
   
   do {
      plvlX--;
      upnnz0.ps--;
      upidx0.ps--;
      upnnz1.ps--;
      upidx1.ps--;
      upnnz2.ps--;
      upidx2.ps--;
      upnnz3.ps--;
      upidx3.ps--;
      lvl = *plvlX;
      
      if( lvl ) {
         register const LSU2X * rnnzLSU3 = (const LSU2X *)&nodes[ *(upnnz3.ps) ];
         register const LSU2X * ridxLSU3 = (const LSU2X *)&nodes[ *(upidx3.ps) ];
         lvl--;
         *plvlX = lvl; plvlX++; 
         *(upidx3.pl) = ridxLSU3->r.l; upidx3.pl++;
         *(upnnz3.pl) = rnnzLSU3->r.l; upnnz3.pl++;
         *(upidx3.pl) = ridxLSU3->l.l; upidx3.pl++;
         *(upnnz3.pl) = rnnzLSU3->l.l; upnnz3.pl++;
         register const LSU2X * rnnzLSU2 = (const LSU2X *)&nodes[ *(upnnz2.ps) ];
         register const LSU2X * ridxLSU2 = (const LSU2X *)&nodes[ *(upidx2.ps) ];
         *plvlX = lvl; plvlX++; 
         *(upidx2.pl) = ridxLSU2->r.l; upidx2.pl++;
         *(upnnz2.pl) = rnnzLSU2->r.l; upnnz2.pl++;
         *(upidx2.pl) = ridxLSU2->l.l; upidx2.pl++;
         *(upnnz2.pl) = rnnzLSU2->l.l; upnnz2.pl++;
                                xnnzLSU1 = (const LSU2X *)&nodes[ *(upnnz1.ps) ];
         register const LSU2X * ridxLSU1 = (const LSU2X *)&nodes[ *(upidx1.ps) ];
         *plvlX = lvl; plvlX++; 
         *(upidx1.pl) = ridxLSU1->r.l; upidx1.pl++;
         *(upnnz1.pl) = xnnzLSU1->r.l; upnnz1.pl++;
         *(upidx1.pl) = ridxLSU1->l.l; upidx1.pl++;
         *(upnnz1.pl) = xnnzLSU1->l.l; upnnz1.pl++;
                                xnnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ];
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         *plvlX = lvl; plvlX++; 
         *(upidx0.pl) = ridxLSU0->r.l; upidx0.pl++;
         *(upnnz0.pl) = xnnzLSU0->r.l; upnnz0.pl++;
         *(upidx0.pl) = ridxLSU0->l.l; upidx0.pl++;
         *(upnnz0.pl) = xnnzLSU0->l.l; upnnz0.pl++;
         
         if( xnnzLSU1->l.s.lo == 0 ) {
            DEBUG_PRINT("... ... in SpMVFlatKernelCOO_HL4_last(), both sides traversal, breaking on empty left tree." );
            break;
         }
      }
      else {
         // leaf node
         register const LSU2X * rnnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ];
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         
         register S32 rowBrk00 = ridxLSU0->l.s.lo;
         register S32 colIdx00 = ridxLSU0->l.s.hi + 1 + colIdx01;
                      colIdx01 = ridxLSU0->r.s.hi + 1 + colIdx00;
                      
         if( rowBrk00 ) {
            *y0 = result0;
            result0 = 0.0;
            y0 += rowBrk00;
         }
         register double xv00 = _x[ colIdx00 ];
         register double xv01 = _x[ colIdx01 ];
            
         register double nzv00 = rnnzLSU0->l.d;
         register double nzv01 = rnnzLSU0->r.d;
            
         result0 += xv00 * nzv00 + xv01 * nzv01;
         //////////////////////////////////////////////////////////////////////////////////////////
         register const LSU2X * rnnzLSU1 = (const LSU2X *)&nodes[ *(upnnz1.ps) ];         
         register const LSU2X * ridxLSU1 = (const LSU2X *)&nodes[ *(upidx1.ps) ];
         
         register S32 rowBrk10 = ridxLSU1->l.s.lo;
         register S32 colIdx10 = ridxLSU1->l.s.hi + 1 + colIdx11;
                      colIdx11 = ridxLSU1->r.s.hi + 1 + colIdx10;
         
         
         if( rowBrk10 ) {
            *y1 = result1;
            result1 = 0.0;
            y1 += rowBrk10;
         }
         register double xv10 = _x[ colIdx10 ];
         register double xv11 = _x[ colIdx11 ];
            
         register double nzv10 = rnnzLSU1->l.d;
         register double nzv11 = rnnzLSU1->r.d;
            
         result1 += xv10 * nzv10 + xv11 * nzv11;
         //////////////////////////////////////////////////////////////////////////////////////////
         register const LSU2X * rnnzLSU2 = (const LSU2X *)&nodes[ *(upnnz2.ps) ];
         register const LSU2X * ridxLSU2 = (const LSU2X *)&nodes[ *(upidx2.ps) ];
         
         register S32 rowBrk20 = ridxLSU2->l.s.lo;
         register S32 colIdx20 = ridxLSU2->l.s.hi + 1 + colIdx21;
                      colIdx21 = ridxLSU2->r.s.hi + 1 + colIdx20;
         
         
         if( rowBrk20 ) {
            *y2 = result2;
            result2 = 0.0;
            y2 += rowBrk20;
         }
         register double xv20 = _x[ colIdx20 ];
         register double xv21 = _x[ colIdx21 ];
            
         register double nzv20 = rnnzLSU2->l.d;
         register double nzv21 = rnnzLSU2->r.d;
            
         result2 += xv20 * nzv20 + xv21 * nzv21;
         //////////////////////////////////////////////////////////////////////////////////////////
         register const LSU2X * rnnzLSU3 = (const LSU2X *)&nodes[ *(upnnz3.ps) ];         
         register const LSU2X * ridxLSU3 = (const LSU2X *)&nodes[ *(upidx3.ps) ];
         
         register S32 rowBrk30 = ridxLSU3->l.s.lo;
         register S32 colIdx30 = ridxLSU3->l.s.hi + 1 + colIdx31;
                      colIdx31 = ridxLSU3->r.s.hi + 1 + colIdx30;
         
         
         if( rowBrk30 ) {
            *y3 = result3;
            result3 = 0.0;
            y3 += rowBrk30;
         }
         register double xv30 = _x[ colIdx30 ];
         register double xv31 = _x[ colIdx31 ];
            
         register double nzv30 = rnnzLSU3->l.d;
         register double nzv31 = rnnzLSU3->r.d;
            
         result3 += xv30 * nzv30 + xv31 * nzv31;
         //////////////////////////////////////////////////////////////////////////////////////////
      }
   } while( plvlX != final );
   
   *y1 = result1;
   *y2 = result2;
   *y3 = result3;
   
   if( xnnzLSU0->r.s.lo )
   while( plvlX != final ) {
      plvlX--;
      upnnz0.ps--;
      upidx0.ps--;
      lvl = *plvlX;
      
      if( lvl ) {
         // internal node
         lvl--;
         *plvlX = lvl; plvlX++;
         *plvlX = lvl; plvlX++;
                                xnnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ];
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         *plvlX = lvl; plvlX++; 
         *plvlX = lvl; plvlX++;
         *(upidx0.pl) = ridxLSU0->r.l; upidx0.pl++;
         *(upnnz0.pl) = xnnzLSU0->r.l; upnnz0.pl++;
         *(upidx0.pl) = ridxLSU0->l.l; upidx0.pl++;
         *(upnnz0.pl) = xnnzLSU0->l.l; upnnz0.pl++;
         if( xnnzLSU0->l.s.lo == 0 ) {
            DEBUG_PRINT("... ... in SpMVFlatKernelCOO_HL4_last(), left side only traversal, breaking on empty left tree." );
            break;
         }
      }
      else {
         // leaf node
         register const LSU2X * rnnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ];
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         
         register S32 rowBrk00 = ridxLSU0->l.s.lo;
         register S32 colIdx00 = ridxLSU0->l.s.hi + 1 + colIdx01;
                      colIdx01 = ridxLSU0->r.s.hi + 1 + colIdx00;
                      
         if( rowBrk00 ) {
            *y0 = result0;
            result0 = 0.0;
            y0 += rowBrk00;
         }
         register double xv00 = _x[ colIdx00 ];
         register double xv01 = _x[ colIdx01 ];
            
         register double nzv00 = rnnzLSU0->l.d;
         register double nzv01 = rnnzLSU0->r.d;
            
         result0 += xv00 * nzv00 + xv01 * nzv01;
      }
   }
   
   *y0 = result0;
   
   DEBUG_PRINT( "... ... leaving SpMVFlatKernelCOO_HL4_last() " );
   #endif
   return 0;
}
inline void * SpMVFlatKernelCOO_HL4_full( void * vpargs ) {
   #ifdef __HL4__
   DEBUG_PRINT( "... ... entering SpMVFlatKernelCOO_HL4_full() " );
   
   COOArgs * pargs = (COOArgs*) vpargs;
   const U64 hgt = pargs->hgt;
   const U64 * _nodes = pargs->_nodes;

   const U32 nnzValsRootNodePtr0 = pargs->nnzValsRootNodePtr[ 0 ];
   const U32 colIdxsRootNodePtr0 = pargs->colIdxsRootNodePtr[ 0 ];
   const U32 nnzValsRootNodePtr1 = pargs->nnzValsRootNodePtr[ 1 ];
   const U32 colIdxsRootNodePtr1 = pargs->colIdxsRootNodePtr[ 1 ];
   const U32 nnzValsRootNodePtr2 = pargs->nnzValsRootNodePtr[ 2 ];
   const U32 colIdxsRootNodePtr2 = pargs->colIdxsRootNodePtr[ 2 ];
   const U32 nnzValsRootNodePtr3 = pargs->nnzValsRootNodePtr[ 3 ];
   const U32 colIdxsRootNodePtr3 = pargs->colIdxsRootNodePtr[ 3 ];
   
   const double * _x = pargs->_x;
   double * _y = pargs->_y;
   double * y0 = _y;
   double * y1 = _y;
   double * y2 = _y;
   double * y3 = _y;
   y0--;
   y1--;
   y2--;
   y3--;
   
   const U32 maxHgt   = 20;
   const U32 numAlloc = 4 * maxHgt;
   U32 lvl = hgt;

   U32 levelsStackX[ numAlloc ];
   U32 nnzPtrStack0[ numAlloc ];
   U32 idxPtrStack0[ numAlloc ];
   U32 nnzPtrStack1[ numAlloc ];
   U32 idxPtrStack1[ numAlloc ];
   U32 nnzPtrStack2[ numAlloc ];
   U32 idxPtrStack2[ numAlloc ];
   U32 nnzPtrStack3[ numAlloc ];
   U32 idxPtrStack3[ numAlloc ];
   
   const S32 * nodes = (const S32 *) _nodes;
   register S32 colIdx01 = 0;
   register S32 colIdx11 = 0;
   register S32 colIdx21 = 0;
   register S32 colIdx31 = 0;
   register double result0 = 0.0;
   register double result1 = 0.0;
   register double result2 = 0.0;
   register double result3 = 0.0;
   
   U32 * plvlX  = &levelsStackX[0];

   PLS upidx0;
   PLS upnnz0;
   PLS upidx1;
   PLS upnnz1;
   PLS upidx2;
   PLS upnnz2;
   PLS upidx3;
   PLS upnnz3;
   
   upidx0.ps = &idxPtrStack0[0];
   upnnz0.ps = &nnzPtrStack0[0];
   upidx1.ps = &idxPtrStack1[0];
   upnnz1.ps = &nnzPtrStack1[0];
   upidx2.ps = &idxPtrStack2[0];
   upnnz2.ps = &nnzPtrStack2[0];
   upidx3.ps = &idxPtrStack3[0];
   upnnz3.ps = &nnzPtrStack3[0];
   
   const U32 * final = plvlX;
   
   const LSU2X * xnnzLSU0 = (const LSU2X *)&nodes[ nnzValsRootNodePtr0 ];
   const LSU2X * xidxLSU0 = (const LSU2X *)&nodes[ colIdxsRootNodePtr0 ];
   const LSU2X * xnnzLSU1 = (const LSU2X *)&nodes[ nnzValsRootNodePtr1 ];
   const LSU2X * xidxLSU1 = (const LSU2X *)&nodes[ colIdxsRootNodePtr1 ];
   const LSU2X * xnnzLSU2 = (const LSU2X *)&nodes[ nnzValsRootNodePtr2 ];
   const LSU2X * xidxLSU2 = (const LSU2X *)&nodes[ colIdxsRootNodePtr2 ];
   const LSU2X * xnnzLSU3 = (const LSU2X *)&nodes[ nnzValsRootNodePtr3 ];
   const LSU2X * xidxLSU3 = (const LSU2X *)&nodes[ colIdxsRootNodePtr3 ];

   *plvlX = lvl; plvlX++;
   *plvlX = lvl; plvlX++;
   *plvlX = lvl; plvlX++;
   *plvlX = lvl; plvlX++;
   
   *(upnnz0.pl) = xnnzLSU0->r.l; upnnz0.pl++;
   *(upnnz0.pl) = xnnzLSU0->l.l; upnnz0.pl++;
   *(upidx0.pl) = xidxLSU0->r.l; upidx0.pl++;
   *(upidx0.pl) = xidxLSU0->l.l; upidx0.pl++;

   *(upnnz1.pl) = xnnzLSU1->r.l; upnnz1.pl++;
   *(upnnz1.pl) = xnnzLSU1->l.l; upnnz1.pl++;
   *(upidx1.pl) = xidxLSU1->r.l; upidx1.pl++;
   *(upidx1.pl) = xidxLSU1->l.l; upidx1.pl++;

   *(upnnz2.pl) = xnnzLSU2->r.l; upnnz2.pl++;
   *(upnnz2.pl) = xnnzLSU2->l.l; upnnz2.pl++;
   *(upidx2.pl) = xidxLSU2->r.l; upidx2.pl++;
   *(upidx2.pl) = xidxLSU2->l.l; upidx2.pl++;

   *(upnnz3.pl) = xnnzLSU3->r.l; upnnz3.pl++;
   *(upnnz3.pl) = xnnzLSU3->l.l; upnnz3.pl++;
   *(upidx3.pl) = xidxLSU3->r.l; upidx3.pl++;
   *(upidx3.pl) = xidxLSU3->l.l; upidx3.pl++;
   
   
   do {
      plvlX--;
      upnnz0.ps--;
      upidx0.ps--;
      upnnz1.ps--;
      upidx1.ps--;
      upnnz2.ps--;
      upidx2.ps--;
      upnnz3.ps--;
      upidx3.ps--;
      lvl = *plvlX;
      
      if( lvl ) {
         // internal node
         register const LSU2X * rnnzLSU3 = (const LSU2X *)&nodes[ *(upnnz3.ps) ];
         register const LSU2X * ridxLSU3 = (const LSU2X *)&nodes[ *(upidx3.ps) ];
         lvl--;
         *plvlX = lvl; plvlX++; 
         *(upidx3.pl) = ridxLSU3->r.l; upidx3.pl++;
         *(upnnz3.pl) = rnnzLSU3->r.l; upnnz3.pl++;
         *(upidx3.pl) = ridxLSU3->l.l; upidx3.pl++;
         *(upnnz3.pl) = rnnzLSU3->l.l; upnnz3.pl++;
         register const LSU2X * rnnzLSU2 = (const LSU2X *)&nodes[ *(upnnz2.ps) ];
         register const LSU2X * ridxLSU2 = (const LSU2X *)&nodes[ *(upidx2.ps) ];
         *plvlX = lvl; plvlX++; 
         *(upidx2.pl) = ridxLSU2->r.l; upidx2.pl++;
         *(upnnz2.pl) = rnnzLSU2->r.l; upnnz2.pl++;
         *(upidx2.pl) = ridxLSU2->l.l; upidx2.pl++;
         *(upnnz2.pl) = rnnzLSU2->l.l; upnnz2.pl++;
         register const LSU2X * rnnzLSU1 = (const LSU2X *)&nodes[ *(upnnz1.ps) ];
         register const LSU2X * ridxLSU1 = (const LSU2X *)&nodes[ *(upidx1.ps) ];
         *plvlX = lvl; plvlX++; 
         *(upidx1.pl) = ridxLSU1->r.l; upidx1.pl++;
         *(upnnz1.pl) = rnnzLSU1->r.l; upnnz1.pl++;
         *(upidx1.pl) = ridxLSU1->l.l; upidx1.pl++;
         *(upnnz1.pl) = rnnzLSU1->l.l; upnnz1.pl++;
         register const LSU2X * rnnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ];
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         *plvlX = lvl; plvlX++; 
         *(upidx0.pl) = ridxLSU0->r.l; upidx0.pl++;
         *(upnnz0.pl) = rnnzLSU0->r.l; upnnz0.pl++;
         *(upidx0.pl) = ridxLSU0->l.l; upidx0.pl++;
         *(upnnz0.pl) = rnnzLSU0->l.l; upnnz0.pl++;
      }
      else {
         // leaf node
         register const LSU2X * rnnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ];
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         
         register S32 rowBrk00 = ridxLSU0->l.s.lo;
         register S32 colIdx00 = ridxLSU0->l.s.hi + 1 + colIdx01;
                      colIdx01 = ridxLSU0->r.s.hi + 1 + colIdx00;
                      
         if( rowBrk00 ) {
            *y0 = result0;
            result0 = 0.0;
            y0 += rowBrk00;
         }
         register double xv00 = _x[ colIdx00 ];
         register double xv01 = _x[ colIdx01 ];
            
         register double nzv00 = rnnzLSU0->l.d;
         register double nzv01 = rnnzLSU0->r.d;
            
         result0 += xv00 * nzv00 + xv01 * nzv01;
         //////////////////////////////////////////////////////////////////////////////////////////
         register const LSU2X * rnnzLSU1 = (const LSU2X *)&nodes[ *(upnnz1.ps) ];         
         register const LSU2X * ridxLSU1 = (const LSU2X *)&nodes[ *(upidx1.ps) ];
         
         register S32 rowBrk10 = ridxLSU1->l.s.lo;
         register S32 colIdx10 = ridxLSU1->l.s.hi + 1 + colIdx11;
                      colIdx11 = ridxLSU1->r.s.hi + 1 + colIdx10;
         
         
         if( rowBrk10 ) {
            *y1 = result1;
            result1 = 0.0;
            y1 += rowBrk10;
         }
         register double xv10 = _x[ colIdx10 ];
         register double xv11 = _x[ colIdx11 ];
            
         register double nzv10 = rnnzLSU1->l.d;
         register double nzv11 = rnnzLSU1->r.d;
            
         result1 += xv10 * nzv10 + xv11 * nzv11;
         //////////////////////////////////////////////////////////////////////////////////////////
         register const LSU2X * rnnzLSU2 = (const LSU2X *)&nodes[ *(upnnz2.ps) ];
         register const LSU2X * ridxLSU2 = (const LSU2X *)&nodes[ *(upidx2.ps) ];
         
         register S32 rowBrk20 = ridxLSU2->l.s.lo;
         register S32 colIdx20 = ridxLSU2->l.s.hi + 1 + colIdx21;
                      colIdx21 = ridxLSU2->r.s.hi + 1 + colIdx20;
         
         
         if( rowBrk20 ) {
            *y2 = result2;
            result2 = 0.0;
            y2 += rowBrk20;
         }
         register double xv20 = _x[ colIdx20 ];
         register double xv21 = _x[ colIdx21 ];
            
         register double nzv20 = rnnzLSU2->l.d;
         register double nzv21 = rnnzLSU2->r.d;
            
         result2 += xv20 * nzv20 + xv21 * nzv21;
         //////////////////////////////////////////////////////////////////////////////////////////
         register const LSU2X * rnnzLSU3 = (const LSU2X *)&nodes[ *(upnnz3.ps) ];         
         register const LSU2X * ridxLSU3 = (const LSU2X *)&nodes[ *(upidx3.ps) ];
         
         register S32 rowBrk30 = ridxLSU3->l.s.lo;
         register S32 colIdx30 = ridxLSU3->l.s.hi + 1 + colIdx31;
                      colIdx31 = ridxLSU3->r.s.hi + 1 + colIdx30;
         
         
         if( rowBrk30 ) {
            *y3 = result3;
            result3 = 0.0;
            y3 += rowBrk30;
         }
         register double xv30 = _x[ colIdx30 ];
         register double xv31 = _x[ colIdx31 ];
            
         register double nzv30 = rnnzLSU3->l.d;
         register double nzv31 = rnnzLSU3->r.d;
            
         result3 += xv30 * nzv30 + xv31 * nzv31;
         //////////////////////////////////////////////////////////////////////////////////////////
      }   
   } while( plvlX != final );
   
   *y0 = result0;
   *y1 = result1;
   *y2 = result2;
   *y3 = result3;
   
   // DEBUG_PRINT( "leaving SpMVFlatKernelCOO_HL4_full() " );
   #endif
   return 0;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
inline void * SpMVFlatKernelCOO_HL2_last( void * vpargs ) {
   #ifdef __HL2__
   DEBUG_PRINT( "... ... entering SpMVFlatKernelCOO_HL2_last() " );
   COOArgs * pargs = (COOArgs*) vpargs;
   const U64 hgt = pargs->hgt;
   const U64 * _nodes = pargs->_nodes;

   const U32 nnzValsRootNodePtr0 = pargs->nnzValsRootNodePtr[ 0 ];
   const U32 colIdxsRootNodePtr0 = pargs->colIdxsRootNodePtr[ 0 ];
   const U32 nnzValsRootNodePtr1 = pargs->nnzValsRootNodePtr[ 1 ];
   const U32 colIdxsRootNodePtr1 = pargs->colIdxsRootNodePtr[ 1 ];
   
   const double * _x = pargs->_x;
   double * _y = pargs->_y;
   double * y0 = _y;
   double * y1 = _y;
   y0--;
   y1--;

   const LSU2X * nnzLSU0;
   const LSU2X * idxLSU0;
   const LSU2X * nnzLSU1;
   const LSU2X * idxLSU1;
   
   const U32 maxHgt   = 20;
   const U32 numAlloc = 4 * maxHgt;

   U32 levelsStackX[ numAlloc ];
   U32 nnzPtrStack0[ numAlloc ];
   U32 idxPtrStack0[ numAlloc ];
   U32 nnzPtrStack1[ numAlloc ];
   U32 idxPtrStack1[ numAlloc ];
   
   const S32 * nodes = (const S32 *) _nodes;
   register S32 colIdx01 = 0;
   register S32 colIdx11 = 0;
   register double result0 = 0.0;
   register double result1 = 0.0;

   PLS upidx0;
   PLS upidx1;
   PLS upnnz0;
   PLS upnnz1;
   
   upidx0.ps = &idxPtrStack0[0];
   upidx1.ps = &idxPtrStack1[0];
   upnnz0.ps = &nnzPtrStack0[0];
   upnnz1.ps = &nnzPtrStack1[0];

   U32 * plvlX  = &levelsStackX[0];

   const U32 * final = plvlX;

   nnzLSU0 = (const LSU2X *)&nodes[ nnzValsRootNodePtr0 ];
   idxLSU0 = (const LSU2X *)&nodes[ colIdxsRootNodePtr0 ];
   nnzLSU1 = (const LSU2X *)&nodes[ nnzValsRootNodePtr1 ];
   idxLSU1 = (const LSU2X *)&nodes[ colIdxsRootNodePtr1 ];

   *plvlX = hgt; plvlX++;
   *plvlX = hgt; plvlX++;
   *plvlX = hgt; plvlX++;
   *plvlX = hgt; plvlX++;
   
   *(upnnz0.pl) = nnzLSU0->r.l; upnnz0.pl++;
   *(upnnz0.pl) = nnzLSU0->l.l; upnnz0.pl++;
   *(upidx0.pl) = idxLSU0->r.l; upidx0.pl++;
   *(upidx0.pl) = idxLSU0->l.l; upidx0.pl++;

   *(upnnz1.pl) = nnzLSU1->r.l; upnnz1.pl++;
   *(upnnz1.pl) = nnzLSU1->l.l; upnnz1.pl++;
   *(upidx1.pl) = idxLSU1->r.l; upidx1.pl++;
   *(upidx1.pl) = idxLSU1->l.l; upidx1.pl++;
   // U32 lvl;
   do {
      plvlX--;
      upnnz0.ps--;
      upidx0.ps--;
      upnnz1.ps--;
      upidx1.ps--;
      register U32 lvl = *plvlX;
      //lvl = *plvlX;
      
      if( lvl ) {
         
         lvl--;
                
         //*plvlX = lvl; plvlX++;
         //*plvlX = lvl; plvlX++;
         //*plvlX = lvl; plvlX++;
         *plvlX = lvl; plvlX++;
         register const LSU2X * ridxLSU1 = (const LSU2X *)&nodes[ *(upidx1.ps) ];
         //*plvlX = lvl; plvlX++;
         *(upidx1.pl) = ridxLSU1->r.l; upidx1.pl++;
         *(upidx1.pl) = ridxLSU1->l.l; upidx1.pl++;
         
         *plvlX = lvl; plvlX++;
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         //*plvlX = lvl; plvlX++;
         *(upidx0.pl) = ridxLSU0->r.l; upidx0.pl++;
         *(upidx0.pl) = ridxLSU0->l.l; upidx0.pl++;

         *plvlX = lvl; plvlX++;
         nnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ]; 
         //*plvlX = lvl; plvlX++;
         *(upnnz0.pl) = nnzLSU0->r.l; upnnz0.pl++;
         *(upnnz0.pl) = nnzLSU0->l.l; upnnz0.pl++;

         *plvlX = lvl; plvlX++;
         register const LSU2X * rnnzLSU1 = (const LSU2X *)&nodes[ *(upnnz1.ps) ]; 
         //*plvlX = lvl; plvlX++;
         *(upnnz1.pl) = rnnzLSU1->r.l; upnnz1.pl++;
         *(upnnz1.pl) = rnnzLSU1->l.l; upnnz1.pl++;         
         if( rnnzLSU1->l.s.lo == 0 ) {
            DEBUG_PRINT("... ... in SpMVFlatKernelCOO_HL2_last(), both sides traversal, breaking on empty left tree." );
            // if( rnnzLSU0->r.s.lo == 0 ) { goto end; }
            break;
         }
      }
      else {
         // leaf node
         register const LSU2X * rnnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ];
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         register const LSU2X * rnnzLSU1 = (const LSU2X *)&nodes[ *(upnnz1.ps) ];
         register const LSU2X * ridxLSU1 = (const LSU2X *)&nodes[ *(upidx1.ps) ];
         
         register S32 rowBrk00 = ridxLSU0->l.s.lo;
         register S32 colIdx00 = ridxLSU0->l.s.hi + 1 + colIdx01;
                      colIdx01 = ridxLSU0->r.s.hi + 1 + colIdx00;
         
         if( rowBrk00 ) {
            *y0 = result0;
            result0 = 0.0;
            y0 += rowBrk00;
         }

         register S32 rowBrk10 = ridxLSU1->l.s.lo;
         register S32 colIdx10 = ridxLSU1->l.s.hi + 1 + colIdx11;
                      colIdx11 = ridxLSU1->r.s.hi + 1 + colIdx10;

         if( rowBrk10 ) {
            *y1 = result1;
            result1 = 0.0;
            y1 += rowBrk10;
         }

         register double xv00 = _x[ colIdx00 ];
         register double xv01 = _x[ colIdx01 ];
         register double xv10 = _x[ colIdx10 ];
         register double xv11 = _x[ colIdx11 ];

         register double nzv00 = rnnzLSU0->l.d;
         register double nzv01 = rnnzLSU0->r.d;
         register double nzv10 = rnnzLSU1->l.d;
         register double nzv11 = rnnzLSU1->r.d;

         result0 += xv00 * nzv00 + xv01 * nzv01;
         result1 += xv10 * nzv10 + xv11 * nzv11;
         // DEBUG_PRINT( "... ... <<< leaf node done >>> " );
      }
   } while( plvlX != final );
   
   if( nnzLSU0->r.s.lo == 0 )
   DEBUG_PRINT("... ... in SpMVFlatKernelCOO_HL2_last(): 1/2 way through." );
   while( plvlX != final ) {
      plvlX--;
      upnnz0.ps--;
      upidx0.ps--;
      register U32 lvl = *plvlX;
      
      if( lvl ) {
         // internal node
         lvl--;
         *plvlX = lvl; plvlX++;
         *plvlX = lvl; plvlX++;
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         *(upidx0.pl) = ridxLSU0->r.l; upidx0.pl++;
         *(upidx0.pl) = ridxLSU0->l.l; upidx0.pl++;
         
         *plvlX = lvl; plvlX++;
         *plvlX = lvl; plvlX++;
         register const LSU2X * rnnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ];
         *(upnnz0.pl) = rnnzLSU0->r.l; upnnz0.pl++;
         *(upnnz0.pl) = rnnzLSU0->l.l; upnnz0.pl++;
         if( rnnzLSU0->l.s.lo == 0 ) {
            DEBUG_PRINT("... ... in SpMVFlatKernelCOO_HL2_last(), left side only traversal, breaking on empty left tree." );
            break;
         }
      }
      else {
         // leaf node
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         
         register S32 rowBrk00 = ridxLSU0->l.s.lo;
         register S32 colIdx00 = ridxLSU0->l.s.hi + 1 + colIdx01;
                      colIdx01 = ridxLSU0->r.s.hi + 1 + colIdx00;
         
         if( rowBrk00 ) {
            *y0 = result0;
            result0 = 0.0;
            y0 += rowBrk00;
         }

         register double xv00 = _x[ colIdx00 ];
         register double xv01 = _x[ colIdx01 ];
         
         register const LSU2X * rnnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ];
         register double nzv00 = rnnzLSU0->l.d;
         register double nzv01 = rnnzLSU0->r.d;

         result0 += xv00 * nzv00 + xv01 * nzv01;
         // DEBUG_PRINT( "... ... <<< leaf node done >>> " );         
      }
   }
   *y1 = result1;
   *y0 = result0;
   
   DEBUG_PRINT( "... ... leaving SpMVFlatKernelCOO_HL2_last() " );
   #endif
   return 0;
}
inline void * SpMVFlatKernelCOO_HL2_full( void * vpargs ) {
   #ifdef __HL2__
   DEBUG_PRINT( "... ... entering SpMVFlatKernelCOO_HL2_full() " );
   COOArgs * pargs = (COOArgs*) vpargs;
   const U64 hgt = pargs->hgt;
   const U64 * _nodes = pargs->_nodes;

   const U32 nnzValsRootNodePtr0 = pargs->nnzValsRootNodePtr[ 0 ];
   const U32 colIdxsRootNodePtr0 = pargs->colIdxsRootNodePtr[ 0 ];
   const U32 nnzValsRootNodePtr1 = pargs->nnzValsRootNodePtr[ 1 ];
   const U32 colIdxsRootNodePtr1 = pargs->colIdxsRootNodePtr[ 1 ];
   
   const double * _x = pargs->_x;
   double * _y = pargs->_y;
   double * y0 = _y;
   double * y1 = _y;
   y0--;
   y1--;
   
   const LSU2X * nnzLSU0;
   const LSU2X * idxLSU0;
   const LSU2X * nnzLSU1;
   const LSU2X * idxLSU1;
   
   const U32 maxHgt   = 20;
   const U32 numAlloc = 4 * maxHgt;
   U32 lvl = hgt;

   U32 levelsStackX[ numAlloc ];
   U32 nnzPtrStack0[ numAlloc ];
   U32 idxPtrStack0[ numAlloc ];
   U32 nnzPtrStack1[ numAlloc ];
   U32 idxPtrStack1[ numAlloc ];
   
   const S32 * nodes = (const S32 *) _nodes;
   register S32 colIdx01 = 0;
   register S32 colIdx11 = 0;
   register double result0 = 0.0;
   register double result1 = 0.0;
   
   U32 * plvlX  = &levelsStackX[0];

   PLS upidx0;
   PLS upidx1;
   PLS upnnz0;
   PLS upnnz1;
   
   upidx0.ps = &idxPtrStack0[0];
   upidx1.ps = &idxPtrStack1[0];
   upnnz0.ps = &nnzPtrStack0[0];
   upnnz1.ps = &nnzPtrStack1[0];

   const U32 * final = plvlX;
   
   nnzLSU0 = (const LSU2X *)&nodes[ nnzValsRootNodePtr0 ];
   idxLSU0 = (const LSU2X *)&nodes[ colIdxsRootNodePtr0 ];
   nnzLSU1 = (const LSU2X *)&nodes[ nnzValsRootNodePtr1 ];
   idxLSU1 = (const LSU2X *)&nodes[ colIdxsRootNodePtr1 ];

   *plvlX = lvl; plvlX++;
   *plvlX = lvl; plvlX++;
   *plvlX = lvl; plvlX++;
   *plvlX = lvl; plvlX++;
   
   *(upnnz0.pl) = nnzLSU0->r.l; upnnz0.pl++;
   *(upnnz0.pl) = nnzLSU0->l.l; upnnz0.pl++;
   *(upidx0.pl) = idxLSU0->r.l; upidx0.pl++;
   *(upidx0.pl) = idxLSU0->l.l; upidx0.pl++;

   *(upnnz1.pl) = nnzLSU1->r.l; upnnz1.pl++;
   *(upnnz1.pl) = nnzLSU1->l.l; upnnz1.pl++;
   *(upidx1.pl) = idxLSU1->r.l; upidx1.pl++;
   *(upidx1.pl) = idxLSU1->l.l; upidx1.pl++;
   
   do {
      plvlX--;
      upnnz0.ps--;
      upidx0.ps--;
      upnnz1.ps--;
      upidx1.ps--;
      lvl = *plvlX;
      
      if( lvl ) {
         // internal node
         register const LSU2X * rnnzLSU1 = (const LSU2X *)&nodes[ *(upnnz1.ps) ]; 
         register const LSU2X * ridxLSU1 = (const LSU2X *)&nodes[ *(upidx1.ps) ];
         
         lvl--;         
         *plvlX = lvl; plvlX++; 
         *(upidx1.pl) = ridxLSU1->r.l; upidx1.pl++;
         *(upnnz1.pl) = rnnzLSU1->r.l; upnnz1.pl++;
         
         
         *plvlX = lvl; plvlX++;
         *(upidx1.pl) = ridxLSU1->l.l; upidx1.pl++;
         *(upnnz1.pl) = rnnzLSU1->l.l; upnnz1.pl++;

         register const LSU2X * rnnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ]; 
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         *plvlX = lvl; plvlX++;
         *(upidx0.pl) = ridxLSU0->r.l; upidx0.pl++;
         *(upnnz0.pl) = rnnzLSU0->r.l; upnnz0.pl++;
         
         *plvlX = lvl; plvlX++;
         *(upidx0.pl) = ridxLSU0->l.l; upidx0.pl++;
         *(upnnz0.pl) = rnnzLSU0->l.l; upnnz0.pl++;
      }
      else {
         // leaf node
         register const LSU2X * rnnzLSU0 = (const LSU2X *)&nodes[ *(upnnz0.ps) ];
         register const LSU2X * ridxLSU0 = (const LSU2X *)&nodes[ *(upidx0.ps) ];
         
         register S32 rowBrk00 = ridxLSU0->l.s.lo;
         register S32 colIdx00 = ridxLSU0->l.s.hi + 1 + colIdx01;
                      colIdx01 = ridxLSU0->r.s.hi + 1 + colIdx00;
                      
         if( rowBrk00 ) {
            *y0 = result0;
            result0 = 0.0;
            y0 += rowBrk00;
         }
         register double xv00 = _x[ colIdx00 ];
         register double xv01 = _x[ colIdx01 ];
            
         register double nzv00 = rnnzLSU0->l.d;
         register double nzv01 = rnnzLSU0->r.d;
            
         result0 += xv00 * nzv00 + xv01 * nzv01;

         register const LSU2X * rnnzLSU1 = (const LSU2X *)&nodes[ *(upnnz1.ps) ];         
         register const LSU2X * ridxLSU1 = (const LSU2X *)&nodes[ *(upidx1.ps) ];
         
         register S32 rowBrk10 = ridxLSU1->l.s.lo;
         register S32 colIdx10 = ridxLSU1->l.s.hi + 1 + colIdx11;
                      colIdx11 = ridxLSU1->r.s.hi + 1 + colIdx10;
         
         
         if( rowBrk10 ) {
            *y1 = result1;
            result1 = 0.0;
            y1 += rowBrk10;
         }
         register double xv10 = _x[ colIdx10 ];
         register double xv11 = _x[ colIdx11 ];
            
         register double nzv10 = rnnzLSU1->l.d;
         register double nzv11 = rnnzLSU1->r.d;
            
         result1 += xv10 * nzv10 + xv11 * nzv11;
      }
   } while( plvlX != final );
   
   *y0 = result0;
   *y1 = result1;
   #endif
   return 0;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
inline void * SpMVFlatKernelCOO_HL1_last( void * vpargs ) {
   #ifdef __HL1__
   DEBUG_PRINT( "... ... entering SpMVFlatKernelCOO() " );
   COOArgs * pargs = (COOArgs*) vpargs;
   const U64 hgt = pargs->hgt;
   const U64 * _nodes = pargs->_nodes;

   const U32 nnzValsRootNodePtr = pargs->nnzValsRootNodePtr[ 0 ];
   const U32 colIdxsRootNodePtr = pargs->colIdxsRootNodePtr[ 0 ];

   const double * _x = pargs->_x;
   double * _y = pargs->_y;
   double * y = _y;
   y--;
   
   const LSU * nnzLSU;
   const LSU * idxLSU;

   
   const U32 maxHgt   = 20;
   const U32 numAlloc = 4 * maxHgt;
   U32 lvl = hgt;

   U32 levelsStack[ numAlloc ];
   U32 nnzPtrStack[ numAlloc ];
   U32 idxPtrStack[ numAlloc ];
   
   const S32 * nodes = (const S32 *) _nodes;
   register S32 colIdx1 = 0;
   register double result = 0.0;
   
   U32 * plvl  = &levelsStack[0];
   U32 * pidx  = &idxPtrStack[0];
   U32 * pnnz  = &nnzPtrStack[0];
   
   nnzLSU = (const LSU *)&nodes[ nnzValsRootNodePtr ]; nnzLSU++; 
   idxLSU = (const LSU *)&nodes[ colIdxsRootNodePtr ]; idxLSU++;

   *plvl = lvl; plvl++; *pidx = idxLSU->s.hi; pidx++; *pnnz = nnzLSU->s.hi; pnnz++;
   *plvl = lvl; plvl++; *pidx = idxLSU->s.lo; pidx++; *pnnz = nnzLSU->s.lo; pnnz++; idxLSU--; nnzLSU--;
   *plvl = lvl; plvl++; *pidx = idxLSU->s.hi; pidx++; *pnnz = nnzLSU->s.hi; pnnz++;
   *plvl = lvl; plvl++; *pidx = idxLSU->s.lo; pidx++; *pnnz = nnzLSU->s.lo; pnnz++;
   
   do {
      plvl--;
      pnnz--;
      pidx--;
      lvl = *plvl;
      
      if( lvl ) {
         // internal node
         nnzLSU = (const LSU *)&nodes[ *pnnz ]; 
         idxLSU = (const LSU *)&nodes[ *pidx ];
         idxLSU++;
         nnzLSU++;
         
         lvl--;
         
         *plvl = lvl; plvl++; *pidx = idxLSU->s.hi; pidx++; *pnnz = nnzLSU->s.hi; pnnz++;
         *plvl = lvl; plvl++; *pidx = idxLSU->s.lo; pidx++; *pnnz = nnzLSU->s.lo; pnnz++; idxLSU--; nnzLSU--;
         *plvl = lvl; plvl++; *pidx = idxLSU->s.hi; pidx++; *pnnz = nnzLSU->s.hi; pnnz++;
         *plvl = lvl; plvl++; *pidx = idxLSU->s.lo; pidx++; *pnnz = nnzLSU->s.lo; pnnz++;
         if( nnzLSU->s.lo == 0 ) {
            DEBUG_PRINT("... ... in SpMVFlatKernelCOO_last(), breaking on empty left tree." );
            break;
         }
      }
      else {
         // leaf node
         nnzLSU = (const LSU *)&nodes[ *pnnz ];
         idxLSU = (const LSU *)&nodes[ *pidx ];
         
         register S32 rowBrk0 = idxLSU->s.lo;
         register S32 colIdx0 = idxLSU->s.hi + 1 + colIdx1; idxLSU++;
                      colIdx1 = idxLSU->s.hi + 1 + colIdx0;
         
         register double xv0 = _x[ colIdx0 ];
         register double xv1 = _x[ colIdx1 ];
         
         if( rowBrk0 ) {
            *y = result;
            result = 0.0;
            y += rowBrk0;
         }
         register double nzv0 = nnzLSU->d; nnzLSU++;
         register double nzv1 = nnzLSU->d;
         
         result += xv0 * nzv0 + xv1 * nzv1;
      }
   } while( true );
   
   *y = result;
   DEBUG_PRINT( "... ... leaving SpMVFlatKernelCOO_last() " );
   #endif
   return 0;
}
inline void * SpMVFlatKernelCOO_HL1_full( void * vpargs ) {
   #ifdef __HL1__
   DEBUG_PRINT( "... ... entering SpMVFlatKernelCOO() " );
   COOArgs * pargs = (COOArgs*) vpargs;
   const U64 hgt = pargs->hgt;
   const U64 * _nodes = pargs->_nodes;

   const U32 nnzValsRootNodePtr = pargs->nnzValsRootNodePtr[ 0 ];
   const U32 colIdxsRootNodePtr = pargs->colIdxsRootNodePtr[ 0 ];

   const double * _x = pargs->_x;
   double * _y = pargs->_y;
   double * y = _y;
   y--;
   
   const LSU * nnzLSU;
   const LSU * idxLSU;

   
   const U32 maxHgt   = 20;
   const U32 numAlloc = 4 * maxHgt;
   U32 lvl = hgt;

   U32 levelsStack[ numAlloc ];
   U32 nnzPtrStack[ numAlloc ];
   U32 idxPtrStack[ numAlloc ];
   
   const S32 * nodes = (const S32 *) _nodes;
   register S32 colIdx1 = 0;
   register double result = 0.0;
   
   U32 * plvl  = &levelsStack[0];
   U32 * pidx  = &idxPtrStack[0];
   U32 * pnnz  = &nnzPtrStack[0];
   const U32 * final = plvl;
   
   nnzLSU = (const LSU *)&nodes[ nnzValsRootNodePtr ]; nnzLSU++; 
   idxLSU = (const LSU *)&nodes[ colIdxsRootNodePtr ]; idxLSU++;

   *plvl = lvl; plvl++; *pidx = idxLSU->s.hi; pidx++; *pnnz = nnzLSU->s.hi; pnnz++;
   *plvl = lvl; plvl++; *pidx = idxLSU->s.lo; pidx++; *pnnz = nnzLSU->s.lo; pnnz++; idxLSU--; nnzLSU--;
   *plvl = lvl; plvl++; *pidx = idxLSU->s.hi; pidx++; *pnnz = nnzLSU->s.hi; pnnz++;
   *plvl = lvl; plvl++; *pidx = idxLSU->s.lo; pidx++; *pnnz = nnzLSU->s.lo; pnnz++;
   
   do {
      plvl--;
      pnnz--;
      pidx--;
      lvl = *plvl;
      
      if( lvl ) {
         // internal node
         nnzLSU = (const LSU *)&nodes[ *pnnz ]; 
         idxLSU = (const LSU *)&nodes[ *pidx ];
         idxLSU++;
         nnzLSU++;
         
         lvl--;
         
         *plvl = lvl; plvl++; *pidx = idxLSU->s.hi; pidx++; *pnnz = nnzLSU->s.hi; pnnz++;
         *plvl = lvl; plvl++; *pidx = idxLSU->s.lo; pidx++; *pnnz = nnzLSU->s.lo; pnnz++; idxLSU--; nnzLSU--;
         *plvl = lvl; plvl++; *pidx = idxLSU->s.hi; pidx++; *pnnz = nnzLSU->s.hi; pnnz++;
         *plvl = lvl; plvl++; *pidx = idxLSU->s.lo; pidx++; *pnnz = nnzLSU->s.lo; pnnz++;
         // -- assumes a full tree -- if( nnzLSU->s.lo == 0 ) {
         // -- assumes a full tree --    break; 
         // -- assumes a full tree -- }
      }
      else {
         // leaf node
         nnzLSU = (const LSU *)&nodes[ *pnnz ];
         idxLSU = (const LSU *)&nodes[ *pidx ];
         
         register S32 rowBrk0 = idxLSU->s.lo;
         register S32 colIdx0 = idxLSU->s.hi + 1 + colIdx1; idxLSU++;
                      colIdx1 = idxLSU->s.hi + 1 + colIdx0;
         
         register double xv0 = _x[ colIdx0 ];
         register double xv1 = _x[ colIdx1 ];
         
         if( rowBrk0 ) {
            *y = result;
            result = 0.0;
            y += rowBrk0;
         }
         register double nzv0 = nnzLSU->d; nnzLSU++;
         register double nzv1 = nnzLSU->d;
         
         result += xv0 * nzv0 + xv1 * nzv1;
      }
   } while( plvl != final );
   
   *y = result;
   DEBUG_PRINT(" in SpMVFlatKernelCOO_full(), final = " << final << ", plvl = " << plvl );
   DEBUG_PRINT( "leaving SpMVFlatKernelCOO_full() " );
   #endif
   return 0;
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void * SpMVFlatKernelCOO_HL1_KernelLauncher( void * vpargs ) { return NULL; }
void * SpMVFlatKernelCOO_HL4_KernelLauncher( void * vpargs ) {
   #ifdef __USE_PTHREADS__

   COOArgs * pargs = ( COOArgs * ) vpargs;
   U32 numKernels = pargs->numKernels;
   const bool hasLastKernel = pargs->hasLastKernel;
   
   DEBUG_PRINT( "... inside thread " << pargs->tidx << ", setting thread affinity." );
   SetThreadAffinity( pargs->tidx );
   
   #ifdef __CPUTIME__
   assert( false );
   #endif
   
   DEBUG_PRINT( "... inside thread " << pargs->tidx << ", numKernels = " << numKernels );
   if( hasLastKernel ) {
      numKernels--;
      for( U32 kidx = 0; kidx < numKernels; kidx++ ) {
         
         DEBUG_PRINT( "... inside thread " << pargs->tidx << ", launching kernel " << kidx << " of " << numKernels );
         SpMVFlatKernelCOO_HL4_full( pargs );
         pargs++;
      }
      SpMVFlatKernelCOO_HL4_last( pargs );
   }
   else {
      for( U32 kidx = 0; kidx < numKernels; kidx++ ) {
         
         DEBUG_PRINT( "... inside thread " << pargs->tidx << ", launching kernel " << kidx << " of " << numKernels );
         SpMVFlatKernelCOO_HL4_full( pargs );
         pargs++;
      }
   }
   #endif
   return 0;
}
void * SpMVFlatKernelCOO_HL2_KernelLauncher( void * vpargs ) {

   #ifdef __USE_PTHREADS__

   COOArgs * pargs = ( COOArgs * ) vpargs;
   U32 numKernels = pargs->numKernels;
   const bool hasLastKernel = pargs->hasLastKernel;
   
   DEBUG_PRINT( "... inside thread " << pargs->tidx << ", setting thread affinity." );
   SetThreadAffinity( pargs->tidx );
   
   #ifdef __CPUTIME__
   assert( false );
   #endif
   
   DEBUG_PRINT( "... inside thread " << pargs->tidx << ", numKernels = " << numKernels );
   if( hasLastKernel ) {
      numKernels--;
      for( U32 kidx = 0; kidx < numKernels; kidx++ ) {
         
         DEBUG_PRINT( "... inside thread " << pargs->tidx << ", launching kernel " << kidx << " of " << numKernels );
         SpMVFlatKernelCOO_HL2_full( pargs );
         pargs++;
      }
      SpMVFlatKernelCOO_HL2_last( pargs );
   }
   else {
      for( U32 kidx = 0; kidx < numKernels; kidx++ ) {
         
         DEBUG_PRINT( "... inside thread " << pargs->tidx << ", launching kernel " << kidx << " of " << numKernels );
         SpMVFlatKernelCOO_HL2_full( pargs );
         pargs++;
      }
   }
   #endif
   return 0;   
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void AssignWorkToThreads( const MtxMetaData & mtxdata, const U32 numThreads, const double * x, double * y, CLKS_T clks, COOArgs * args ) {
   
   DEBUG_PRINT( "AssignWorkToThreads()" );

   const bool HL1 = numParts ==   4 ? true : false ;
   const bool HL2 = numParts ==  64 ? true : false ;
   const bool HL4 = numParts == 256 ? true : false ;

   const U32 partsPerKernel = HL1 ? 1 : HL2 ? 2 : 4;

   const S32 nnzLogical     = (( 1 << mtxdata.nnzHgt ) << mtxdata.nnzHgt ) << 1;
   const S32 nnzPerPart     = nnzLogical / numParts;
   const S32 nnzPerKernel   = nnzPerPart * partsPerKernel;
   const U32 partsPerThread = ( numParts + numThreads - 1 ) / numThreads;
   
   S32 nnzRemain  = mtxdata.nnz;
   U32 pidx = 0;
   U32 numAssigned = 0;
   vector< U32 > kernelsPerThread;
   
   for( U32 tidx = 0; tidx < numThreads; tidx++ ) { kernelsPerThread.push_back( 0 ); }
   
   while( nnzRemain > 0 ) {
      for( U32 tidx = 0; tidx < numThreads; tidx++ ) {
         
         const U32 baseIdx = tidx * partsPerThread;
         const U32 argIdx  = baseIdx + numAssigned;
         
         DEBUG_PRINT( "AssignWorkToThreads(): assigning work, nnzRemain = " << SETDEC( 10 ) << nnzRemain << ", pidx = " << SETDEC( 2 ) << pidx << ", tidx = " << tidx << ", argIdx = " << argIdx );
         
         args[ baseIdx ].numKernels++;
         kernelsPerThread[ tidx ]++;
         
         args[ argIdx ].tidx = tidx;
         #ifdef __CPUTIME__
         args[ argIdx ].clks = clks[ tidx ];
         #endif
         args[ argIdx ].nnz  = mtxdata.nnz;
         args[ argIdx ].m    = mtxdata.m;
         args[ argIdx ]._nodes = const_cast<U64*>( mtxdata.idxRoots.nodes );
         
         if( HL1 ) {
            args[ argIdx ].hgt = mtxdata.nnzHgt - 3;
            args[ argIdx ].nnzValsRootNodePtr[ 0 ] = mtxdata.nnzRoots.level1[ pidx ];
            args[ argIdx ].colIdxsRootNodePtr[ 0 ] = mtxdata.idxRoots.level1[ pidx ];
         }
         if( HL2 ) {
            args[ argIdx ].hgt = mtxdata.nnzHgt - 4;
            args[ argIdx ].nnzValsRootNodePtr[ 0 ] = mtxdata.nnzRoots.level2[ pidx + 0 ];
            args[ argIdx ].nnzValsRootNodePtr[ 1 ] = mtxdata.nnzRoots.level2[ pidx + 1 ];
            args[ argIdx ].colIdxsRootNodePtr[ 0 ] = mtxdata.idxRoots.level2[ pidx + 0 ];
            args[ argIdx ].colIdxsRootNodePtr[ 1 ] = mtxdata.idxRoots.level2[ pidx + 1 ];
         }
         if( HL4 ) {
            args[ argIdx ].hgt = mtxdata.nnzHgt - 5;
            args[ argIdx ].nnzValsRootNodePtr[ 0 ] = mtxdata.nnzRoots.level3[ pidx + 0 ];
            args[ argIdx ].nnzValsRootNodePtr[ 1 ] = mtxdata.nnzRoots.level3[ pidx + 1 ];
            args[ argIdx ].nnzValsRootNodePtr[ 2 ] = mtxdata.nnzRoots.level3[ pidx + 2 ];
            args[ argIdx ].nnzValsRootNodePtr[ 3 ] = mtxdata.nnzRoots.level3[ pidx + 3 ];
            args[ argIdx ].colIdxsRootNodePtr[ 0 ] = mtxdata.idxRoots.level3[ pidx + 0 ];
            args[ argIdx ].colIdxsRootNodePtr[ 1 ] = mtxdata.idxRoots.level3[ pidx + 1 ];         
            args[ argIdx ].colIdxsRootNodePtr[ 2 ] = mtxdata.idxRoots.level3[ pidx + 2 ];
            args[ argIdx ].colIdxsRootNodePtr[ 3 ] = mtxdata.idxRoots.level3[ pidx + 3 ];         
         }
         
         args[ argIdx ]._x = x;
         args[ argIdx ]._y = y;

         pidx += partsPerKernel;
         nnzRemain -= nnzPerKernel;
         
         if( nnzRemain <= 0 ) {
            args[ baseIdx ].hasLastKernel = true;
            break;
         }
      }
      numAssigned++;
   }
   
   for( U32 tidx = 0; tidx < numThreads; tidx++ ) {
      DEBUG_PRINT( "thread: " << SETDEC( 2 ) << tidx << ", kernel launches: " << SETDEC( 2 ) << kernelsPerThread[ tidx ] );
   }
   DEBUG_PRINT( "AssignWorkToThreads(): number of partitions assigned = " << pidx );
}
void SpMVCOOThreaded( const MtxMetaData & mtxdata, const double * x, double * y, const U32 numRuns, const U32 numThreads, CLKS_T clks ) {
   
   #ifdef __USE_PTHREADS__

   const bool HL1 = numParts ==   4 ? true : false ;
   const bool HL2 = numParts ==  64 ? true : false ;
   const bool HL4 = numParts == 256 ? true : false ;

   // number of mtx partitions per thread: ceiling( numPartitions / numThreads )
   const U32 partsPerThread   = ( numParts + numThreads - 1 ) / numThreads;
   const U32 klaunchPerThread = partsPerThread / numPartsPerKernel;
   const U32 numArgsStructs   = numThreads * partsPerThread;
   
   pthread_t * threads = (pthread_t*) calloc( numThreads,     sizeof( pthread_t ) );
   COOArgs   * args    = (COOArgs*)   calloc( numArgsStructs, sizeof( COOArgs   ) );
   
   DEBUG_PRINT( "SpMVCOOThreaded(): calling AssignWorkToThreads()" );
   AssignWorkToThreads( mtxdata, numThreads, x, y, clks, &args[0] );

   COOArgs * pargs = &args[ 0 ];
   U32 threadsJoined = 0;
   
   void* (*klaunch)( void* );
   if     ( HL1 ) { klaunch = SpMVFlatKernelCOO_HL1_KernelLauncher; }
   else if( HL2 ) { klaunch = SpMVFlatKernelCOO_HL2_KernelLauncher; }
   else if( HL4 ) { klaunch = SpMVFlatKernelCOO_HL4_KernelLauncher; }
   
   DEBUG_PRINT( "" );
   cout << "numThreads: " << numThreads << endl;
   cout << "mtxPartitionsPerThread:  " << partsPerThread   << endl;
   cout << "kernelLaunchesPerThread: " << klaunchPerThread << endl;
   
   for( U32 tidx = 0; tidx < numThreads; tidx++ ) {

      DEBUG_PRINT( " <<< launching thread " << tidx << " >>> " );
      pthread_create( &threads[ tidx ], NULL, klaunch, pargs );

      pargs += partsPerThread;
      #ifdef __SERIALIZE_PTHREADS__
      pthread_join( threads[ tidx ], NULL );
      threadsJoined++;
      #endif
   }
   for( U32 tidx = threadsJoined; tidx < numThreads; tidx++ ) {
      DEBUG_PRINT( "joining thread " << tidx );
      pthread_join( threads[ tidx ], NULL );
   }

   #ifdef __CHECK_SPMV_RESULT__
   PrintPTRVectorToMMFile( x, mtxdata.n, false, "x_hcoo.mtx" );
   PrintPTRVectorToMMFile( y, mtxdata.m, false, "y_hcoo.mtx" );
   #endif

   #ifdef __SPMV__
   // begin warm cache runs.
   for( U32 runIdx = 0; runIdx < numRuns; runIdx++ ) {

      COOArgs * pargs = &args[ 0 ];
      threadsJoined = 0;
      
      #ifndef __CPUTIME__
      clock_gettime( CLOCK_MONOTONIC, clks );
      clks++;
      #endif
      for( U32 tidx = 0; tidx < numThreads; tidx++ ) {

         DEBUG_PRINT( " <<< launching thread " << tidx << " >>> " );
         pthread_create( &threads[ tidx ], NULL, klaunch, pargs );

         pargs += partsPerThread;
         #ifdef __SERIALIZE_PTHREADS__
         pthread_join( threads[ tidx ], NULL );
         threadsJoined++;
         #endif
      }
      for( U32 tidx = threadsJoined; tidx < numThreads; tidx++ ) {
         DEBUG_PRINT( "joining thread " << tidx );
         pthread_join( threads[ tidx ], NULL );
      }
      #ifndef __CPUTIME__
      clock_gettime( CLOCK_MONOTONIC, clks );
      clks++;
      #endif
   }
   free( threads );
   free( args );
   #endif // __SPMV__
   #endif // __USE_PTHREADS__
}
void SpMVCOOSerial( const MtxMetaData & mtxdata, const DAGLevelPointers & nnzRoots, const DAGLevelPointers & idxRoots, const double * x, double * y, const U32 numRuns, timespec ** clks ) {
   
   U32 i;
   
   COOArgs args[2];
   for( i=0; i < 2; i++ ) {
      args[i].nnz = 0;
      args[i].m   = mtxdata.m;
      args[i].hgt = mtxdata.nnzHgt - 2;
      args[i]._nodes = const_cast<U64*>(nnzRoots.nodes);
      args[i].nnzValsRootNodePtr[0] = nnzRoots.level0[ 2*i + 0 ];
      args[i].colIdxsRootNodePtr[0] = idxRoots.level0[ 2*i + 0 ];
      args[i].nnzValsRootNodePtr[1] = nnzRoots.level0[ 2*i + 1 ];
      args[i].colIdxsRootNodePtr[1] = idxRoots.level0[ 2*i + 1 ];
      args[i]._x = x;
      args[i]._y = y;
   }
   SpMVFlatKernelCOO_HL2_last( &args[0] );
   SpMVFlatKernelCOO_HL2_last( &args[1] );

   #ifdef __CHECK_SPMV_RESULT__
   PrintPTRVectorToMMFile( x, mtxdata.n, false, "x_hcoo.mtx" );
   PrintPTRVectorToMMFile( y, mtxdata.m, false, "y_hcoo.mtx" );
   #endif

   #ifdef __SPMV__
   // begin warm cache runs.
   for( U32 j = 0; j < numRuns; j++ ) {
      clock_gettime( CLOCK_THREAD_CPUTIME_ID, clks[0] ); clks[0]++;
      SpMVFlatKernelCOO_HL2_last( &args[0] );
      SpMVFlatKernelCOO_HL2_last( &args[1] );
      clock_gettime( CLOCK_THREAD_CPUTIME_ID, clks[0] ); clks[0]++;
   }
   #endif
}
void SpMVCOO( MtxMetaData & mtxdata, const U64 * _nodes, const double * x, double * y, const U32 numRuns, const U32 numThreads, CLKS_T clks ) {

   U32 i;
   U32 j;
   
   const U8 offMap[4] = { 1, 0, 3, 2 };
   const U64 mask = 0xFFFFFFFF;
   const U32 nnzValsRootNodePtr = ( mask & _nodes[3] ) << 2;
   const U32 colIdxsRootNodePtr = ( mask & _nodes[5] ) << 2;
   
   const S32 * nodes = (S32*)&_nodes[0];
      
   mtxdata.nnzRoots.nodes = const_cast<U64*>( _nodes );
   mtxdata.idxRoots.nodes = const_cast<U64*>( _nodes );
   
   for( i=0; i <  4; i++ ) {
      mtxdata.nnzRoots.level0[ i ] = nodes[ nnzValsRootNodePtr + offMap[ i ] ];
      mtxdata.idxRoots.level0[ i ] = nodes[ colIdxsRootNodePtr + offMap[ i ] ];
   }
   for( i=0; i <  4; i++ ) {
      for( j=0; j < 4; j++ ) {
         U32 idx = i*4 + j;
         mtxdata.nnzRoots.level1[ idx ] = nodes[ mtxdata.nnzRoots.level0[ i ] + offMap[ j ] ];
         mtxdata.idxRoots.level1[ idx ] = nodes[ mtxdata.idxRoots.level0[ i ] + offMap[ j ] ];
      }
   }
   for( i=0; i < 16; i++ ) {
      for( j=0; j < 4; j++ ) {
         U32 idx = i*4 + j;
         mtxdata.nnzRoots.level2[ idx ] = nodes[ mtxdata.nnzRoots.level1[ i ] + offMap[ j ] ];
         mtxdata.idxRoots.level2[ idx ] = nodes[ mtxdata.idxRoots.level1[ i ] + offMap[ j ] ];
      }
   }
   for( i=0; i < 64; i++ ) {
      for( j=0; j < 4; j++ ) {
         U32 idx = i*4 + j;
         mtxdata.nnzRoots.level3[ idx ] = nodes[ mtxdata.nnzRoots.level2[ i ] + offMap[ j ] ];
         mtxdata.idxRoots.level3[ idx ] = nodes[ mtxdata.idxRoots.level2[ i ] + offMap[ j ] ];
      }      
   }
   assert( mtxdata.nnzHgt == mtxdata.idxHgt );

   #ifdef __USE_PTHREADS__
   SpMVCOOThreaded( mtxdata, x, y, numRuns, numThreads, clks );
   #else
   assert( false );
   // SpMVCOOSerial( mtxdata, x, y, numRuns, clks );
   #endif
   
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void allocclks( const U32 numRuns, const U32 numThreads, vector< timespec * > & clks ) {
   // [ n_threads ][ n_runs ]
   
   for( U32 tidx = 0; tidx < numThreads; tidx++ ) {

      const U32 numClksPerThread  = 2;
      const U32 numClksForAllRuns = numRuns * numClksPerThread;

      timespec * clk_alloc = (timespec *) calloc( numClksForAllRuns, sizeof( timespec ) );

      ASSERT( clk_alloc != NULL );

      clks.push_back( clk_alloc );
   }
}
void freeclks( vector< timespec * > & clks ) {
   
   for( U32 i = 0; i < clks.size(); i++ ) {
      free( clks[i] );
   }
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void SpMVWrapperQTS( const U32 numRuns, const U32 numThreads, const MatrixProperties & mtxProps, U64 * nodes ) {

   const U64   m = mtxProps.m;
   const U64   n = mtxProps.n;
   
   double * y = (double*) MALLOC( m * sizeof(double) );
   double * x = (double*) MALLOC( n * sizeof(double) );

   timespec * clks = (timespec*) calloc( 2 * numRuns, sizeof( timespec ) );
   
   // fill x vector with random values
   // initialize y vector to zeros
   #ifdef __CHECK_SPMV_RESULT__
   srand( 0 );
   #endif
   for( U32 i=0; i < n; i++ ) { x[i] = GenRandDouble(); }
   for( U32 i=0; i < m; i++ ) { y[i] = 0.0;             }
   
   DEBUG_PRINT_VEC( x, MIN( 10, n), "x" );
   DEBUG_PRINT( "Calling SpMVQTS()" );
   DEBUG_PRINT( "" );
   
   SpMVQTS( mtxProps, nodes, x, y, numRuns, numThreads, clks );
   
   DEBUG_PRINT( "" );
   DEBUG_PRINT_VEC( y, MIN( 10, m), "y" );
   
   #ifdef __SPMV__
   AnalyzeRunTimeWallTime( clks, numRuns );
   #endif
   free( clks );
   FREE( (char*) x );
   FREE( (char*) y );
}
void SpMVWrapperCSR( const U32 numRuns, const U32 numThreads, const MatrixProperties & mtxProps, U64 * nodes ) {

   const U64   m = mtxProps.m;
   const U64   n = mtxProps.n;
   const U64 nnz = mtxProps.nnz;
   const U64 nnzHgt = mtxProps.dagHgts[ 0 ] - 1;
   const U64 idxHgt = mtxProps.dagHgts[ 1 ] - 1;
   const U64 ptrHgt = mtxProps.dagHgts[ 2 ] - 1;
   
   double * y = (double*) MALLOC( m * sizeof(double) );
   double * x = (double*) MALLOC( n * sizeof(double) );

   vector< timespec * > clks;
   allocclks( numRuns, numThreads, clks );
      
   // fill x vector with random values
   // initialize y vector to zeros
   #ifdef __CHECK_SPMV_RESULT__
   srand( 0 );
   #endif
   for( U32 i=0; i < n; i++ ) { x[i] = GenRandDouble(); }
   for( U32 i=0; i < m; i++ ) { y[i] = 0.0;             }
   
   DEBUG_PRINT_VEC( x, MIN( 10, n), "x" );
   DEBUG_PRINT( "Calling SpMVCSR( " << nnzHgt << ", " << idxHgt << ", " << ptrHgt << ", " << m << ", " << n << ", " << nnz << ", nodes, x, y, " << numRuns << ", clks )" );
   DEBUG_PRINT( "" );
   
   SpMVCSR( nnzHgt, idxHgt, ptrHgt, m, n, nnz, nodes, x, y, numRuns, clks[0] );
   
   DEBUG_PRINT( "" );
   DEBUG_PRINT_VEC( y, MIN( 10, m), "y" );
   
   #ifdef __SPMV__
   AnalyzeRunTimeCPUTime( clks, numRuns, numThreads );
   #endif
   freeclks( clks );
   FREE( (char*) x );
   FREE( (char*) y );
}
void SpMVWrapperCOO( const U32 numRuns, const U32 numThreads, const MatrixProperties & mtxProps, U64 * nodes ) {

   const U64   m = mtxProps.m;
   const U64   n = mtxProps.n;
   const U64 nnz = mtxProps.nnz;
   const U64 nnzHgt = mtxProps.dagHgts[ 0 ] - 1;
   const U64 idxHgt = mtxProps.dagHgts[ 1 ] - 1;
   
   MtxMetaData mtxdata;
   mtxdata.m = m;
   mtxdata.n = n;
   mtxdata.nnzHgt = nnzHgt;
   mtxdata.idxHgt = idxHgt;
   mtxdata.nnz = nnz;
   
   // cache-line aligned malloc:
   double * y = (double*) MALLOC( ( m + 2 ) * sizeof( double ) );
   double * x = (double*) MALLOC( ( n + 0 ) * sizeof( double ) );
         
   #ifdef __CPUTIME__
   vector< timespec * > clks;
   allocclks( numRuns, numThreads, clks );
   #else
   timespec * clks = (timespec*) calloc( 2 * numRuns, sizeof( timespec ) );
   #endif
      
   // fill x vector with random values
   // initialize y vector to zeros
   #ifdef __CHECK_SPMV_RESULT__
   srand( 0 );
   #endif
   for( U32 i=0; i < ( n + 0 ); i++ ) { x[i] = GenRandDouble(); }
   for( U32 i=0; i < ( m + 2 ); i++ ) { y[i] = 0.0;             }
   
   DEBUG_PRINT_VEC( x, MIN( 10, n), "x" );
   DEBUG_PRINT( "Calling SpMVCOO( " << nnzHgt << ", " << idxHgt << ", " << m << ", " << n << ", " << nnz << ", nodes, x, y, " << numRuns << ", clks )" );
   DEBUG_PRINT( "" );
   
   SpMVCOO( mtxdata, nodes, x, y+1, numRuns, numThreads, clks );
   
   DEBUG_PRINT( "" );
   DEBUG_PRINT_VEC( y, MIN( 10, m), "y" );
   
   #ifdef __SPMV__
   #ifdef __CPUTIME__
   AnalyzeRunTimeCPUTime( clks, numRuns, numThreads );
   freeclks( clks );
   #else
   AnalyzeRunTimeWallTime( clks, numRuns );
   free( clks );
   #endif // __CPUTIME__
   #endif // __SPMV__
   
   FREE( (char*) x );
   FREE( (char*) y );
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void COOSequence( const std::string & mtxFileName, const U32 numRuns, const U32 numThreads ) {
   MatrixProperties mtxProps;
   
   DEBUG_PRINT( "Calling LoadDagMtx( " << mtxFileName << ", 4, mtxProps );" );
   U64 * nodes = LoadDagMtx( mtxFileName, 4, mtxProps );
   DEBUG_PRINT( "Done loading dag." );
   
   DEBUG_PRINT( "Calling SpMVWrapper( " << numRuns << ", mtxProps, nodes );" );
   DEBUG_PRINT( "" );
   SpMVWrapperCOO( numRuns, numThreads, mtxProps, nodes );
         
   #ifdef __HUGE_PAGES__
   struct shmid_ds unused_buffer;
   printf( "Deleting shared segment -- return code %i\n", shmctl( shmidNodes, IPC_RMID, &unused_buffer ) );
   HugeFree( shmidNodes );
   #else
   DEBUG_PRINT( "Calling FREE( (char*)nodes );" );
   FREE( (char*)nodes );
   #endif
}
void CSRSequence( const std::string & mtxFileName, const U32 numRuns, const U32 numThreads ) {
   MatrixProperties mtxProps;
   
   DEBUG_PRINT( "Calling LoadDagMtxCSR( " << mtxFileName << ", mtxProps );" );
   U64 * nodes = LoadDagMtx( mtxFileName, 4, mtxProps );
   DEBUG_PRINT( "Done loading dag." );
   
   DEBUG_PRINT( "Calling SpMVWrapper( " << numRuns << ", mtxProps, nodes );" );
   SpMVWrapperCSR( numRuns, numThreads, mtxProps, nodes );
         
   DEBUG_PRINT( "Calling FREE( (char*)nodes );" );
   FREE( (char*)nodes );
}
void QTSSequence( const std::string & mtxFileName, const U32 numRuns, const U32 numThreads ) {
   MatrixProperties mtxProps;
   
   DEBUG_PRINT( "Calling LoadDagMtx( " << mtxFileName << ", 4, mtxProps );" );
   U64 * nodes = LoadDagMtx( mtxFileName, 4, mtxProps );
   DEBUG_PRINT( "Done loading dag." );
   
   DEBUG_PRINT( "Calling SpMVWrapper( " << numRuns << ", mtxProps, nodes );" );
   SpMVWrapperQTS( numRuns, numThreads, mtxProps, nodes );

   DEBUG_PRINT( "Calling FREE( (char*)nodes );" );
   FREE( (char*)nodes );
}
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
void ParseCmdLineArgs( int argc, char ** argv, U32 & numRuns, U32 & numThreads ) {
   if( argc < 2 ) {
      // useful message
      printf( "usage: %s <qtsFile> [numRuns] [numThreads]\n", argv[0] );
      exit( -1 );
   }
   else if( argc == 3 ) {
      // only numRuns specified; keep default value of 1 thread
      numRuns = atoi( argv[2] );
      if( numRuns < 1 ) {
         printf( "must enter numRuns > 0\n" );
         exit( -1 );
      }
   }
   else if( argc > 3) {
      // numRuns & numThreads specified
      numRuns    = atoi( argv[2] );
      numThreads = atoi( argv[3] );
      
      const U32 numParts = 64;
      const U32 numPartsPerKernel = 2;
      
      assert( numParts % numPartsPerKernel == 0 );
      const U32 maxNumThreads = numParts / numPartsPerKernel;
      
      if( numRuns < 1 ) {
         printf( "must enter numRuns > 0\n" );
         exit( -1 );
      }
      if( numThreads < 1 ) {
         printf( "must enter numThreads > 0\n" );
         exit( -1 );
      }
      if( numThreads > maxNumThreads ) {
         printf( "maximum number of threads = %d\n", maxNumThreads );
         exit( -1 );
      }
   }
}
int main( int argc, char ** argv ) {

   U32 numRuns    = 1;
   U32 numThreads = 1;

   // SetProcessSocketAffinity( 0 ); <== handled by numactl --membind=[0|1]
   ParseCmdLineArgs( argc, argv, numRuns, numThreads );
   // SetSchedulerToFT_FIFO();
   
   DEBUG_PRINT( "" );
   #ifdef __USE_PTHREADS__
   #ifdef __SPMV__
   printf( "=== dagSpMV threaded impl ===\n" );
   #endif
   #else
   #ifdef __SPMV__
   printf( "=== dagSpMV standard impl ===\n" );
   #endif
   #endif
   
   std::string mtxFileName = argv[1];
      
   #ifdef __QTS__
   QTSSequence( mtxFileName, numRuns, numThreads );
   #endif
   
   #ifdef __CSR__
   CSRSequence( mtxFileName, numRuns, numThreads );
   #endif

   #ifdef __COO__
   COOSequence( mtxFileName, numRuns, numThreads );
   #endif

   DEBUG_PRINT( "Exiting main()" );
   DEBUG_PRINT( "" );
   return 0;
}
